#ifndef RSIM_VECTOR_H
#define RSIM_VECTOR_H

/* --------------------------------------------------------------------------*
* File: Vector.h                                                             *
* Authors: Kishor Bhalerao                                                   *
* Email : kishor8dm@gmail.com                                                *
* Contributors:                                                              *
* Email:                                                                     *
* copyright (c) 2010 Authors.                                                *
*                                                                            *
* This program is free software: you can redistribute it and/or modify it    *
* under the terms of the GNU General Public License as published by the Free *
* Software Foundation, either version 3 of the License, or any later version.*
*                                                                            *
* This program is distributed in the hope that it will be useful, but        *
* WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY *
* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License   *
* for more details. (http://www.gnu.org/licenses/)                           *
* -------------------------------------------------------------------------- */

#include "CData.h"
#include "Utils.h"

#include <cstdlib>

namespace RSIM{

template<class T>class Matrix_;

template<class T>
class Vector_:public cData<T>{
        protected:
                bool isColVec_;
        public:
                /** Creates a NULL Vector_. Ensure that space is allocated
                by calling Vector_::setLength() before using the variable.*/
                explicit Vector_(const bool& isColVec=true);
                
                /** Deep copy */
                Vector_(const Vector_<T>& V);
                
                /** Deep copy */
                explicit Vector_(const Matrix_<T>& M);
                
                /** 
                @param Length size of the vector
                @param isColVec column vector is default
                */
                explicit Vector_(const int& Length,const bool& isColVec=true);
                
                /**
                @param Length Length of vector
                @param val value of each element
                @param isColVec Column vector is default 
                */
                Vector_(const int& Length, const T& val, const bool& isColVec=true);
                
                /** Careful. No bound checks on the data */
                Vector_(const int& Length, const T *dataPtr,const bool& isColVec=true);
                
                /** @return length of the vector */
                int dim()const;
                
                /** @return i^th element of the vector. 0 based indexing */
                const T& getElement(const int& i)const;
                
                /** Set all values of this vector to zero. If vector is empty, this does nothing.*/
                void Zeros();
                
                /**
                @param ptr Array from which data is to be copied.
                @param index Vector[index] = ptr[0], Vector[index+1] = ptr[1] and so on.
                @param num Number of values to be copied
                */
                void copy(const T *ptr, const int& index, const int& num);
                
                /** Use this to allocate space for an empty vector */
                void resize(const int& NumElements,const bool& isColVec=true);
                
                /** set i^th element to value v */
                void setElement(const int& i,const T& v);
                
                /** toggles the isColVec variable. */
                void Transpose();
                
                /** true if number of elements = 0 */
                bool isNull()const;
                
                /** True if column vector */ 
                bool& isColVector();
                
                /** True if column vector */
                bool isColVector()const;
                
                Vector_<T> operator+(const Vector_<T>& V)const;
                
                /** Returns a column vector*/
                Vector_<T> extractNElements(const int& n)const;
                
                void operator+=(const Vector_<T>& V);
                
                Vector_<T> operator-(const Vector_<T>& V)const;
                
                Vector_<T> operator+(const Matrix_<T>& M)const;
                
                Vector_<T> operator-(const Matrix_<T>& M)const;
                
                /** vector x scalar */
                Vector_<T> operator*(const T& val)const;
                
                /** vector x matrix */
                Matrix_<T> operator*(const Matrix_<T>& M)const;
                
                /** vector x vector */
                Matrix_<T> operator*(const Vector_<T>& V)const;
                
                /** Transpose of a Matrix_ */
                Vector_<T> operator~(void) const;	
                
                /** @return writable reference */
                T& operator () (const int& i);
                
                /** @return constant reference */
                T& operator () (const int& i) const;
                
                /** @return writable reference */
                T& operator [] (const int& i);
                
                /** @return constant reference */
                const T& operator [] (const int& i)const;
                
                /** Deep copy */
                void operator=(const Vector_<T>& V);
                
                /** Deep copy */
                void operator=(const Matrix_<T>& M);
                
                template <class U>
                friend std::ostream& operator<<(ostream& func, const Vector_<U>& V);
                
                template <class U>
                friend Vector_<U> operator*(const double& val, const Vector_<U>& V);			
                
                template <class U>
                friend std::istream& operator>>(istream& func, RSIM::Vector_<U>& V);

}; // class Vector_

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T> Vector_<T>::extractNElements(const int& n)const{
        Vector_<T> v(n);
        double *ptr = v.wPtr();
        for(int i=0;i<n;ptr[i]=this->data_[i++]);
        return v;
}

////////////////////////////////////////////////////////////////////////

template<class T>
inline int Vector_<T>::dim()const{return this->nele_;}

////////////////////////////////////////////////////////////////////////

template<class T>
inline const T& Vector_<T>::getElement(const int& i)const{return this->data_[i];}

////////////////////////////////////////////////////////////////////////

template<class T>
inline void Vector_<T>::setElement(const int& i,const T& v){this->data_[i]=v;}

////////////////////////////////////////////////////////////////////////

template<class T>
inline void Vector_<T>::Transpose(){isColVec_= !isColVec_;}

////////////////////////////////////////////////////////////////////////

template<class T>
inline const T& Vector_<T>::operator[](const int& i)const{return this->data_[i];}

////////////////////////////////////////////////////////////////////////

template<class T>
inline T& Vector_<T>::operator[](const int& i){return this->data_[i];}

////////////////////////////////////////////////////////////////////////

template<class T>
inline T& Vector_<T>::operator()(const int& i){return this->data_[i];}

////////////////////////////////////////////////////////////////////////

template<class T>
inline T& Vector_<T>::operator()(const int& i)const{return this->data_[i];}

////////////////////////////////////////////////////////////////////////

template<class T>
inline bool Vector_<T>::isColVector()const{return isColVec_;}

////////////////////////////////////////////////////////////////////////

template<class T>
inline bool& Vector_<T>::isColVector(){return isColVec_;}

////////////////////////////////////////////////////////////////////////

template<class T>
inline bool Vector_<T>::isNull()const{if(this->nele_)return false; else return true;}

////////////////////////////////////////////////////////////////////////

template<class T>
inline Vector_<T> Vector_<T>::operator~(void) const {
        Vector_<T> op = *this;
        op.isColVector() = !isColVec_;
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T>
inline void Vector_<T>::resize(const int& NumElements,const bool& isColVec){
        this->setLength(NumElements);
        isColVec_=isColVec;
}

template <class T>
inline void Vector_<T>::Zeros(){
        for(int i=0;i<this->nele_;this->data_[i++]=0);
}

////////////////////////////////////////////////////////////////////////

template<class T>
inline Vector_<T>::Vector_(const bool& isColVec):cData<T>(),isColVec_(isColVec){
}

////////////////////////////////////////////////////////////////////////

template<class T>
inline Vector_<T>::Vector_(const Vector_<T>& V):cData<T>(V),isColVec_(V.isColVector()){
}

////////////////////////////////////////////////////////////////////////

template<class T>
inline Vector_<T>::Vector_(const int& Length, const T *dataPtr,const bool& isColVec):
cData<T>(Length,dataPtr),isColVec_(isColVec){}

////////////////////////////////////////////////////////////////////////

template<class T>
inline Vector_<T>::Vector_(const int& Length,const bool& isColVec):cData<T>(Length),isColVec_(isColVec){
}

////////////////////////////////////////////////////////////////////////

template<class T>
inline Vector_<T>::Vector_(const int& Length, const T& val, const bool& isColVec):cData<T>(Length,val),isColVec_(isColVec){
}

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T> Vector_<T>::operator+(const Vector_<T>& V)const{
        int nele2 = V.getLength();	
        
        Utils::rsim_assert( (this->nele_==nele2) || (isColVec_==V.isColVector()),
                        "Vector_ Vector_::operator+(const Vector_& V)const: Inconsistent Dimensions in addition. exiting...\n"
                        );
                        
        Vector_<T> op(this->nele_);
        const T *v = V.Ptr();
        T *opvec = op.wPtr();
        for(int i=0;i<this->nele_;++i)
                opvec[i]=this->data_[i]+v[i];
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T>
void Vector_<T>::operator+=(const Vector_<T>& V){	
        #ifdef RSIM_SAFE
        int nele2 = V.getLength();
        if((this->nele_!=nele2)||(isColVec_!=V.isColVector())){
                cout<<"\nVector_ Vector_::operator+=(const Vector_& V)const:\n\
                Inconsistent Dimensions in addition. exiting...\n";
                exit(0);
        }
        #endif	
        const T *v = V.Ptr();
        for(int i=0;i<this->nele_;++i)
                this->data_[i] += v[i];
}

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T> Vector_<T>::operator-(const Vector_<T>& V)const{
        int nele2 = V.getLength();
        #ifdef RSIM_SAFE
        if((this->nele_!=nele2)|(isColVec_!=V.isColVector())){
                cout<<"\nVector_ Vector_::operator-(const Vector_& V)const:\n\
                Inconsistent Dimensions in Matrix_<T> addition. exiting...\n";
                exit(0);
        }
        #endif
        Vector_<T> op(this->nele_);
        const T *v = V.Ptr();
        T *opvec = op.wPtr();
        for(int i=0;i<this->nele_;++i)
                opvec[i]=this->data_[i]-v[i];
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T> Vector_<T>::operator+(const Matrix_<T>& M)const{
        int row = M.getNRow();
        int col = M.getNCol();

        if(isColVec_ && (col==1)){
                Vector_<T> op(this->nele_);
                T *opvec = op.wPtr();
                
                const T *m = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        opvec[i]=this->data_[i]+m[i];
                return op;
        }
        else if((!isColVec_) && (row == 1)){
                Vector_<T> op(this->nele_);
                op.isColVector() = false;
                T *opvec = op.wPtr();
                const T *m = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        opvec[i] = this->data_[i]+m[i];
                return op;
        }
        cout<<"\nVector_ Vector_::operator+(const Matrix_<T>& M)const\n\
        Inconsistent Dimensions in addition. exiting...\n";
        exit(0);
}

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T> Vector_<T>::operator-(const Matrix_<T>& M)const{
        int row = M.getNRow();
        int col = M.getNCol();
        if(isColVec_ && (col==1)){
                Vector_<T> op(this->nele_);
                T *opvec = op.wPtr();
                
                const T *m = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        opvec[i]=this->data_[i]-m[i];
                return op;
        }
        else if(!isColVec_ && (row == 1)){
                Vector_<T> op(this->nele_);
                op.isColVector() = false;
                T *opvec = op.wPtr();
                const T *m = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        opvec[i] = this->data_[i]-m[i];
                return op;
        }
        cout<<"\nVector_ Vector_::operator-(const Matrix_<T>& M)const\n\
        Inconsistent Dimensions in addition. exiting...\n";
        exit(0);
}

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T> Vector_<T>::operator*(const T& val)const{
        if(!this->nele_){
                cout<<"\nVector_ Vector_::operator*(const T& val)const\n\
                Attempt to use zero vector. About to exit...\n";
                exit(0);
        }
        Vector_<T> op(this->nele_);
        op.isColVector()=isColVec_;
        
        T *opvec = op.wPtr();
        for(int i=0;i<this->nele_;i++){
                opvec[i]=this->data_[i]*val;
        }
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T>
Matrix_<T> Vector_<T>::operator*(const Matrix_<T>& M)const{
        int Arows, Acols, Brows, Bcols;
        if (isColVec_){
                Arows = this->nele_;
                Acols = 1;
        }
        else{
                Arows = 1;
                Acols = this->nele_;
        }
        Brows = M.getNRow();
        Bcols = M.getNCol();
        #ifdef RSIM_SAFE
        if(Acols!=Brows){
                cout<<"\nMatrix_<T> Vector_::operator*(const Matrix_<T>& M)const\n\
                Incompatible dimensions. About to exit...\n";
                exit(0);
        }
        #endif
        Matrix_<T> op(Arows,Bcols);
        T *opmat = op.wPtr();
        const T *m = M.Ptr();
        
        T sum;
        for(int i = 0;i<Arows;i++){
                for(int j=0;j<Bcols;j++){
                        sum=0.0;
                        for(int k=0;k<Acols;k++)
                                sum = sum+this->data_[i*Acols+k]*m[k*Bcols+j];
                        opmat[i*Bcols+j]=sum;
                }
        }
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T>
Matrix_<T> Vector_<T>::operator*(const Vector_<T>& V)const{
        int Arows, Acols, Brows, Bcols;
        if (isColVec_){
                Arows = this->nele_;
                Acols = 1;
        }
        else{
                Arows = 1;
                Acols = this->nele_;
        }
        if (V.isColVector()){
                Brows = V.getLength();
                Bcols = 1;
        }
        else{
                Brows = 1;
                Bcols = V.getLength();

        }	
        #ifdef RSIM_SAFE
        if(Acols!=Brows){
                cout<<"\nMatrix_<T> Vector_::operator*(const Vector_& V)const\n\
                        Incompatible Matrix_<T> dimensions. About to exit...\n";
                exit(0);
        }
        #endif
        Matrix_<T> op(Arows,Bcols);
        T *opmat = op.wPtr();
        const T *v2 = V.Ptr();
        T sum;
        for(int i = 0;i<Arows;i++){
                for(int j=0;j<Bcols;j++){
                        sum=0.0;
                        for(int k=0;k<Acols;k++)
                                sum = sum+this->data_[i+k]*v2[k+j];
                        opmat[i*Bcols+j]=sum;
                }
        }
        return op;
}

////////////////////////////////////////////////////////////////////////

template<class T>
void Vector_<T>::operator=(const Vector_<T>& V){
        int len = V.getLength();
        const T *ptr = V.Ptr();
        if((this->nele_==len) && (isColVec_==V.isColVector())){
                for(int i=0;i<len;++i)
                        this->data_[i] = ptr[i];
        }
        else if(!this->nele_){
                this->setLength(len);
                isColVec_ = V.isColVector();
                for(int i=0;i<len;++i)
                        this->data_[i] = ptr[i];
        }
        else{
                cout<<"\nvoid Vector_::operator=(const Vector_& V)\n\
                Incompatible vector dimensions. Existing...\n";
                exit(0);
        }
}

////////////////////////////////////////////////////////////////////////

template<class T>
void Vector_<T>::operator=(const Matrix_<T>& M){
        int r = M.getNRow();
        int c = M.getNCol();
        if( !((r-1)*(c-1)) && !(this->nele_-r*c)){
                if(!(r-1))
                        isColVec_ = false;
                else
                        isColVec_ = true;
                const T *ptr = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        this->data_[i] = ptr[i];
        }
        else if( !((r-1)*(c-1)) && !this->nele_){
                if(!(r-1))
                        isColVec_ = false;
                else
                        isColVec_ = true;
                this->setLength(r*c);
                const T *ptr = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        this->data_[i] = ptr[i];
        }
        else{
                cout<<"\nvoid Vector_::operator=(const Matrix_<T>& M)\n\
                Incompatible dimensions, exiting... \n";
                exit(0);
        }
}

////////////////////////////////////////////////////////////////////////

template<class T>
Vector_<T>::Vector_(const Matrix_<T>& M):cData<T>(),isColVec_(true){
        int r = M.getNRow();
        int c = M.getNCol();
        if( !((r-1)*(c-1)) ){
                if(!(r-1))
                        isColVec_ = false;
                this->setLength(r*c);
                const T *ptr = M.Ptr();
                for(int i=0;i<this->nele_;++i)
                        this->data_[i] = ptr[i];
        }
        else{
                cout<<"\nVector_::Vector_(const Matrix_<T>& M)\n\
                either row or col of the matrix M should be 1. Exiting... \n";
                exit(0);
        }
}
////////////////////////////////////////////////////////////////////////

template<class T>
std::ostream& operator<<(ostream& func, const Vector_<T>& V){ 
        if(V.isColVector())
                func<<"~[";
        else		
                func<<"[";
        for(int i=0;i<V.getLength();i++){
                if(i)
                        func<<", ";
                func<<V.getElement(i);
        }
        func<<"]";
        return func;
}

////////////////////////////////////////////////////////////////////////

template <class T>
Vector_<T> operator*(const double& val, const Vector_<T>& V){
        return (V*val);
}

////////////////////////////////////////////////////////////////////////

template <class T>
std::istream& operator>>(istream& func, Vector_<T>& V){
        for(int i=0;i<V.getLength();i++){
                func>>V(i); 
        }
        return func; 
}

////////////////////////////////////////////////////////////////////////

template <class T>
void Vector_<T>::copy(const T *ptr, const int& index, const int& num){
        for(int i=0;i<num;++i)
                this->data_[index+i] = ptr[i];
}

////////////////////////////////////////////////////////////////////////

}// namespace RSIM

#endif
